Method for determining functional volumes for determining biokinetics

ABSTRACT

A method for determining functional volumes for the search for kinetics representing the trend of concentration of a radioactive tracer in an area of biological tissue, being applied to spatial components comprises the following steps applied iteratively according to a Markov chain Monte-Carlo scheme: generation of a set K 0   λ  made up of a set of candidate kinetics associated with values of probability of appearance, depending on the concentration λ of radioactive tracer; a labeling step during which, for each spatial component of index K, the probabilities of selection of the kinetics of the set K 0   λ  are weighted by introducing a function λ k  representative of the concentration of radioactive tracer in this component to obtain a set of indicative values, an indicative value D k  designating the kinetic with which the spatial component K is associated; construction of functional volumes, a functional volume VF j  made up of the set of spatial components which share the same indicative value D k .

The invention relates to a method for determining functional volumes forpositron emission tomography and applies notably to the fields ofpositron emission tomography and pharmacokinetics.

Positron emission tomography, usually referred to by the acronym PET, isa nuclear medicine imaging technology which uses molecules marked with apositron emitting radionucleide to image, in vivo, the molecularinteractions of biological processes with a minimum of perturbations.PET is the modality of reference for molecular imaging because of itsmolecular sensitivity of the order of the femto-mole and of thequantitative character of the measurements that it makes possible. PETimaging is used in particular to determine the functional volumeshereinafter in the description designated by the acronym VF. Afunctional volume VF corresponds to an area consisting of tissuesexhibiting identical metabolic properties.

The state of the art in the field of estimating functional volumesproceeds indirectly by reconstructing, in a first step, a temporalseries of discretized images then by estimating the functional volumesfrom this series of images. Now, reducing the dose of radioactivetracer, also called radiopharmaceutical, contributes to degrading thesignal-to-noise ratio of each of the images to the detriment, in thisindirect approach, of the determination of the functional volumes.

PET is also often used in clinical procedures in a static, 2D or 3Dspatial context, in which the time variable is disregarded, for want oftools allowing for an analysis of the dynamic data that is robust,automatic and validated. This use is significantly sub-optimal. Ineffect, PET molecular imaging is by nature dynamic inasmuch as thedistribution of the radioactive tracer injected into the organismevolves over time as a function of the molecular interactions betweenthe tracer and the tissues. The objective of the interpretation of theexamination can then be summarized as a quantitative and objectivecharacterization of the spatial extensions exhibiting similar temporalpharmacokinetic behaviors, hence the notion of functional volumes. It isquantitative inasmuch as it has to be pronounced on the uncertaintyassociated with the molecular concentration contained in each functionalvolume in terms of quantity of radioactive tracer. It is objectiveinasmuch as the number, the positioning and the extension of thefunctional volumes do not rely on the supervision of an expert havingfirst defined one or more regions of interest.

The state of the art in the dynamic study of biokinetics by PET imagingis based on a temporal segmentation of the information in order toobtain a 3D series. This means that the acquisition is subdivided intotime slices reconstructed independently of one another. Although itallows for a dynamic tracking of the distribution of the radioactivetracer from the series of reconstructed images, this approach presentsthe drawback that only a proportion of the observations is used toreconstruct an image. This statistical loss fatally induces adeterioration of the quality of the reconstruction which is all thegreater when the time slice is short.

Most of the time, the determination of the spatial extension of thefunctional volume is defined a priori, from the anatomical imaging.Then, from the temporal series of images, the concentration of meanactivity inside the functional volume, which in reality corresponds toan anatomical volume, is measured, and measured for each time slice.Existing methods make it possible to determine the extension of thefunctional volume directly from the dynamic images. Source separatingtechniques are also used to determine the distinct molecular processesfrom the reconstructed images. However, the number of kinetics is fixedin these approaches.

A. Reader et al. take up this principle in their paper entitled Jointestimation of dynamic PET images and temporal basis functions usingfully 4D ML-EM, Physics in Medicine and Biology, 51(21):5455-5474, 2006,by using a dictionary of temporal functions to seek the breakdown in afunctional base which best adjusts the dataset. From a limited number oftemporal basis functions, typically about ten or so, it is possible apriori to describe, by their combination, an infinite number of timeactivity curves, usually referred to by the acronym TAC. These 4Dreconstruction methods make it possible to enhance the signal-to-noiseratio of the reconstructed images. All these techniques seek to performa classification, also referred to by the term “clustering”, ofspace-time data of the disintegrations, but proceed indirectly to thisend since a step of reconstructions of voxelized images over the entirespace is necessary. In effect, these techniques require thediscretization of all the spatial and temporal dimensions, resulting inthe manipulation and storage of large volumes of data and considerablecomputation times. Also, the latter approach does not make it possibleto quantify the estimation error without approximation.

Works on reconstruction in emission tomography without recourse to thediscretization of the space have been conducted in a non-parametricBayesian framework. This approach has been extended to the space-timereconstruction as presented in the paper by D. Fall, E. Barat, C.Comtat, T. Dautremer, T. Montagu, and S. Stute entitled Continuousspace-time reconstruction in 4D PET, IEEE Medical Imaging Conference,2011. This paper is called Fall, Barat et al. hereinafter in thedescription. This method addresses the problem of the quantitativedetermination of the reconstruction error of the functional volumes andthe objective of reducing the dose injected into the patient. It furtherproposes a space-time “clustering” in a continuous space directly fromthe observations. However, in this approach, the TACs are all consideredas normalized kinetics. This modeling, although it allows for theautomatic detection of any number of kinetics, is not optimal inasmuchas the discrimination between kinetics uses only their forms and doesnot take into account the concept of concentration, that is to say theamplitude of the TACs expressed in terms of Becquerel per second, thisunit being denoted Bq.s⁻¹. The limitation of this modeling comes fromthe model retained for the clustering of the spatial components with aparticular kinetic. In effect, the latter implements a hierarchicalmodel in which the random distribution of the kinetics, represented byK₀, is itself distributed according to a Pólya tree Dirichlet process.Then, the distribution of the space-time components denoted H is chosenas a Dirichlet process whose central measurement involves K₀. Noinformation concerning the concentration is present in this stochasticmodel.

One aim of the invention is notably to mitigate the abovementioneddrawbacks.

To this end, the subject of the invention is a method for determiningfunctional volumes for the search for kinetics representing the trend ofthe concentration of a radioactive tracer in an area of biologicaltissue. The method is applied to spatial components and comprises thefollowing steps applied iteratively according to a Markov chainMonte-Carlo scheme:

-   -   a step of generation of a set K₀ ^(λ) made up of a set of        candidate kinetics associated with values of probability of        appearance of these kinetics, these values depending on the        concentration λ of radioactive tracer;    -   a labeling step during which, for each spatial component of        index K, the probabilities of selection of the kinetics of the        set K₀ ^(λ) are weighted by introducing therein a function λ_(k)        representative of the concentration of radioactive tracer in        this component so as to obtain a set of indicative values, an        indicative value D_(k) designating the kinetic with which the        spatial component K is associated;    -   a step of construction of functional volumes, a functional        volume VF_(j) being made up of the set of the spatial components        which share the same indicative value D_(k).

According to one aspect of the invention, a step computes averages ofthe functional volumes VF_(j) obtained during the iterative process.

In one embodiment, a step determines statistical estimators concerningthe functional volumes VF_(j) obtained during the iterative process soas to quantify the estimation uncertainty. The variance of thefunctional volumes VF_(j) obtained during the iterative process can thenbe determined. A credible interval can also be determined.

By way of example, the set K₀ ^(λ) is obtained by the determination oftwo sets of scalar variables V_(j) and Γ_(j), these variables relatingto kinetic components of index j, in which V_(j) represents a fixedweighting coefficient independent of the concentration λ and Γ_(j)represents a reference concentration value associated with the kineticof index j.

According to one aspect of the invention, V_(j) is determined by usingthe expression:

$\mspace{20mu} {{V_{j}D},{T\overset{ind}{\sim}{{Beta}\left( {{1 + m_{j}},{v + {\text{?}m_{l}\text{?}\text{indicates text missing or illegible when filed}}}} \right.}}}$

an expression in which:

-   -   for any j, m_(j) designates the number of components of the        distribution H allocated to the j^(th) component of K₀ ^(Λ), H        being a random measurement distributed according to a partially        dependent Dirichlet process, K₀ ^(Λ) being an innumerable        collection of random probability measurements which accepts, as        distribution, a local Dirichlet process;    -   J_(k) _(n) designates the greatest index of the components of K₀        ^(Λ) having at least one component of H allocated;    -   Beta(,) represents the beta probability law.

In one embodiment, Γ_(j) is determined by using the expression:

${\Gamma_{j}V_{j}},D,w,{\overset{\sim}{Z}\overset{ind}{\sim}{U_{j}^{*}\left( \Gamma_{j} \right)}}$with${U_{j}^{\bigstar}\left( \Gamma_{j} \right)} \propto {\left( {a_{j,\Gamma}^{\prime} < \Gamma_{j} < b_{j,\Gamma}^{\prime}} \right) \times {\prod\limits_{k \in _{j}^{+}}^{\;}\; \left( {1 - {V_{j}\left( {{d\left( {\lambda_{k},\Gamma_{j}} \right)} < \psi} \right)}} \right)}}$

expressions in which:

-   -   Γ_(j) designates the parameter characteristic of the        concentration associated with the kinetic j;    -   D designates the set of the variables D_(k) which indicate with        which kinetic of K₀ ^(λ) the spatial component K is associated;    -   D_(j)={k≦k_(n):D_(k)=j} represents the set of the indices K of        the spatial components associated with the kinetic j;    -   D_(j) ⁺={k≦k_(n):D_(k)>j} represents the set of the indices K of        the spatial components associated with a kinetic of index        greater than j;    -   w designates the set of the weightings w_(k) of the spatial        components k in the Dirichlet mixture H;    -   {tilde over (Z)} designates the set of the parameters, average        and covariance matrix,    -   {tilde over (Z)}_(k) associated with the spatial component K of        the Dirichlet mixture H;    -   here indicates an independent random generation for each        parameter Γ_(j);    -   ∞ here signifies “proportional to”;    -   ( ) represents the “gate” function such that        (condition)=1 if condition is true, and 0 otherwise;    -   V_(j) designates a coefficient associated with each kinetic        involved in the weighting of said kinetic in K₀ ^(λ);    -   d( ) represents the Euclidian distance in        : d(x₁,x₂)=x₁ x₂| in which x₁, x₂∈        ;    -   ψ is a threshold chosen a priori which makes it possible to        parameterize K₀ ^(λ).

According to one aspect of the invention, the labeling step isperformed:

-   -   by independently generating for any k≦k_(n):

$\mspace{20mu} {\left( {{D_{k}V},v,Q^{\bigstar},\Gamma,w,\overset{\sim}{Z},C,T} \right)\overset{ind}{\sim}{\sum\limits_{j = 1}^{J^{\bigstar}}\; {p_{j,k}{\delta_{j}( \cdot )}}}}$  with  p_(j, k) ∝ (p_(j)(λ_(k)) > υ_(k))max (p_(j)(λ_(k)), ?)?f_(Q_(j)^(★))(T_(i))  and$\mspace{20mu} {{\sum\limits_{j = 1}^{J^{\bigstar}}\; p_{j,k}} = 1}$?indicates text missing or illegible when filed

-   -   by independently generating for any K, such that k_(n)<k≦k*:

$\mspace{20mu} {\left( {{D_{k}V},v,Q^{\bigstar},\Gamma,w,\overset{\sim}{Z}} \right)\overset{ind}{\sim}{\sum\limits_{j = 1}^{J^{\bigstar}}\; {p_{j,k}{\delta_{j}( \cdot )}}}}$  with   p_(j, k) ∝ (p_(j)(λ_(k)) > υ_(k))max (p_(j)(λ_(k)), ?)  and$\mspace{20mu} {{\sum\limits_{j = 1}^{J^{\bigstar}}\; p_{j,k}} = 1.}$?indicates text missing or illegible when filed

in these expressions,

-   -   k_(n) designates the number of spatial components with which at        least one coincidence event is associated;    -   k* designates the total number of spatial components of the        Dirichlet mixture H for a given iteration;    -   V designates the set of the coefficients V_(j) involved in the        weightings of the kinetics mixture K₀ ^(λ);    -   v represents the set of the auxiliary variables v_(k), used in        the representation of the Dirichlet kinetics process;    -   Q* designates the set of the Pólya trees Q*_(j) characterizing        the kinetics;    -   Γ designates the set of the parameters Γ_(j) characteristic of        the concentration associated with the kinetic j;    -   C designates the set of the classification variables C_(i) which        associate each coincidence event i with a spatial component of        the Dirichlet mixture H;    -   T represents the set of the times of occurrences T_(i) of the        coincidence events;    -   p_(j,k) designates the probability of associating the component        K of the Dirichlet mixture H with the kinetic j;    -   p_(j)(λ) represents the weight function taken at λ associated        with the probability of the kinetic j of K₀ ^(λ);    -   δ_(j)( ) represents the Dirac measurement function such that        δ_(j)(j′)=1 if j=j′ and is 0 otherwise;    -   here indicates an independent random generation for each        classification variable D_(k);    -   max( ) represents the “maximum of” function;    -   ζ designates is a threshold parameter chosen arbitrarily and        involved in the implementation of the “slice sampling” algorithm        of the local Dirichlet kinetics process;    -   {i: C_(i)=k} designates the set of the coincidence events i        which are associated with the spatial component K;    -   f_(Q*) _(j) represents the kinetic j, i.e. the density of the        Pólya tree j of the local Dirichlet mixture of kinetics;    -   J* designates the total number of components of the local        Dirichlet mixture of kinetics for a given iteration.

According to another aspect of the invention, the concentration λ_(k) isobtained by the application of a function φ_(λ) defined such thatλ_(k)=φ_(λ)({tilde over (Z)}_(k), S) with:

$\mspace{20mu} {\lambda_{k} = {\text{?}\left( {\sum\limits_{m = 1}^{\infty}\; {\text{?}{\left( {\text{?}\text{?}} \right)}}} \right)}}$  with      ? ∼ (??)?indicates text missing or illegible when filed

-   -   in which:    -   _({tilde over (x)}) _(k) designates the mathematical expectation        relative to the Gaussian probability law of {tilde over        (x)}_(k), N({tilde over (x)}_(k)|{tilde over (Z)}_(k));    -   N( ) represents the Gaussian probability law.

Also the subject of the invention is a positron emission tomographydevice implementing the method described previously.

The method according to the invention offers the decisive advantage ofincreasing the distances between functional volumes and therebyincreasing the separation quality.

Also, the method according to the invention lies in a novel approachthat makes it possible to shrewdly establish the a posteriori law of theconcentration of each kinetic, thus conferring on the method thepossibility of providing an estimation of the uncertainty on theasymptotic metabolization parameters of a tissue.

Advantageously, the method according to the invention makes it possibleto substantially reduce the dose injected into the patient for afour-dimensional PET imaging while providing functional information to amolecular scale.

A repeated use of PET imaging during therapeutic monitoring can then beenvisaged without any dosimetric prejudice to the patient. A wider usein pediatrics and in translational research can also be considered.

Furthermore, a reduction of the dose of radioactive tracer injectedautomatically brings about a reduction of the exposure to ionizingradiations of the medical personnel involved in preparing and carryingout the examination and of the entourage of the patient.

Other features and advantages of the invention will become apparent fromthe following description given by way of illustration and as anonlimiting example, in light of the attached drawings in which:

FIG. 1 schematically illustrates the various steps of the method fordetermining functional volumes;

FIG. 2 illustrates how the weights and auxiliary variables of thedistribution H of the space-time components are generated;

FIG. 3 schematically illustrates how the spatial parameters used by themethod are generated;

FIG. 4 schematically illustrates one way of determining the local andauxiliary variables of the set K₀ ^(λ);

FIG. 5 illustrates how the kinetics are allocated to the components ofthe set K₀ ^(λ);

FIG. 6 shows how the places of emission are determined and how thecomponents of H are determined;

FIG. 7 represents one way of computing the concentration λ_(k) ofradioactive tracer;

FIG. 1 schematically illustrates the various steps of the method fordetermining functional volumes.

The aim of the method is to determine functional volumes for the searchfor kinetics. It is applied to measurements presented in the form ofcoincidence events 106 and comprises a series of steps appliediteratively 111 according to a Markov chain Monte-Carlo scheme.

In a first step 108, a set K₀ ^(λ) is generated made up of a set ofcandidate kinetics associated with values of probability of appearanceof these kinetics, these values depending on the concentration λ ofradioactive tracer.

In a so-called labeling second step 109, for each spatial component ofindex k 106, the probabilities of selection of the kinetics of the setK_(o) _(λ) are weighted 102 by introducing therein a function λ_(k).This function is representative of the concentration of radioactivetracer of this component so as to obtain a set of indicative values 103.An indicative value D_(k) designating the kinetic with which the spatialcomponent k is associated is determined.

In a third step 110, the functional volumes are constructed, afunctional volume VF_(j), being made up of the set of the spatialcomponents which share the same indicative value D_(k).

The method for determining functional volumes described hereinbelowmakes it possible to mitigate the independence of the model used byFall, Barat et al. with regard to the concentration information. Toachieve this, the concentration of radioactive tracer is estimated atany point of space by marginalization, that is to say by integrationover time of the distribution of space-time activity. This estimatedconcentration is then considered as a random co-variable endogenous tothe problem. The term endogenous underscores the fact that thisco-variable is not external because the latter is determined bybeginning the reconstruction. Provided with this co-variable denoted λ,the proposed method implements a modeling in the distribution of thekinetics dependent on the concentration of radioactive tracer. Thisamounts to replacing the distribution K₀ with a collection ofdistributions K₀ ^(Λ)={K₀ ^(λ): λ∈Λ}. This collection of distributionsthen accepts an independent Dirichlet process as distribution.

In one embodiment, a local Pólya tree Dirichlet process IDP is chosen.One consequence is that the distribution H of the space-time componentsis no longer generated according to a Dirichlet process but according toa distribution referred to hereinbelow in the description as partiallydependent Dirichlet process and designated by the acronym pdDP.

According to one aspect of the invention, the construction of the pdDPis parameterized by a function φ_(λ) which makes it possible to link theconcentration co-variable λ to the parameters of the marginal spatialdistribution obtained by integration over time of the space-timedistribution.

The general problem of reconstruction of a probability distribution G(·)for which only samples of a projection thereof denoted F(·) can beobserved is formulated in the non-parametric Bayesian framework for theinverse problems by using the following expressions:

G˜

F(·, t)=∫_(x) P(·|x)G(dx, t)

Y _(i) , T _(i) |F

F, for i=1, . . . , n   (1)

in which:

-   -   n is the number of samples observed;    -   G(·) is a random probability distribution distributed according        to        , defined over (X×T, σ(X)        σ(T));    -   the distribution G(·) is estimated from observations distributed        according to F(·): (Y, T)′={(Y₁, T₁), . . . , (Y_(n), T_(n))};    -   P(·|x) is the assumed known probability distribution, called        projector, indexed by x and defined over (y, σ(y));    -   the notation σ(ε) corresponds to the sigma-algebra over ε.

The object that is intended to be reconstructed is the spatialdistribution of activity in the field of view of the tomography. In PET,each coincidence event is assigned to a virtual line LOR, LOR being anacronym for “line-of-response”, joining the centers of the twodetectors. A response line LOR of index l is characterized by theparameters defining its position in the space of the detectors. Acorrelation is constructed below between the parameters of the LORs andthe index l such that the variable y_(i) corresponds to the index of theLOR having recorded the i^(th) observed coincidence. More specifically,the observation y_(i) designates the coordinates of the LOR which joinsthe two detectors having recorded the coinciding photons deriving fromthe annihilation of the positrons having taken place at x_(i).

In this context, the terms of the equation (1) correspond to thefollowing interpretation:

-   -   X⊂        ³ corresponds to the region of the reconstruction space;    -   T(        ⁺ corresponds to the acquisition time band;    -   Y_(i) corresponds to the index of the observed LOR;    -   T_(i) corresponds to the time of occurrence of the i^(th) event        observed;    -   P(y=m|x) for m=1, . . . , M is the distribution corresponding to        the projector, that is to say to the probability of recording an        event on the m^(th) detection LOR for a localized emission at x,        that is to say the place of annihilation for any x∈X,

${{\sum\limits_{m = 1}^{M}\; {\left( {y = {mx}} \right)}} = {{\int_{y}^{\;}{\left( {{y}x} \right)}} = 1}}\ $

-   -   G(·) is the space-time distribution of the emissions        corresponding to the recorded events. The set of observations is        truncated because of the limited axial field of the imager.

In order to take account of these truncations in the reconstruction, anLOR of zero index is defined, the latter corresponding to thenon-observable emissions such that y*=y∪{0}. It is then possible todefine the probability of projection P*(y−m|x) over (y*, σ(y*)) and

${\left( {yx} \right)} = {\frac{^{*}\left( {yx} \right)}{\int_{y}^{\;}{^{*}\left( \ {{y}x} \right)}} = \frac{^{*}\left( {yx} \right)}{\sum\limits_{m = 1}^{M}{^{*}\left( {{y - m}x} \right)}}}$

The space-time density of the places of emissions can then be written:

${G^{\bigstar}\left( {x,t} \right)} \propto \frac{G\left( {x,t} \right)}{\int_{y}^{\;}{^{*}\left( \ {{y}x} \right)}}$

an expression in which ∫_(y)P*(dy|x) designates the probability ofrecording a localized emission at x by taking into account the geometryof the system and the physical properties of the detection system.

The expression (2) formally defines the proposed model which makes itpossible to take into account the spatial concentration as endogenousco-variable of the reconstruction problem:

Y _(i) |X _(i)

P(Y _(i) |X _(i))

X _(i) , T _(i) |Z _(i) , Q _(i)

N(X _(i) |Z _(i))×Q _(i)(T _(i))

Z _(i), Q_(i) |H

H

H˜pdDP(θ, G ₀ , K ₀ ⁻, φ_(x))

K₀ ^(Λ) ={K ₀ ^(λ): λ∈Λ}˜IDP(v, PT _(L)(A, Q ₀), γ₀, ψ)   (2)

In this expression,

-   -   H is a random measurement distributed according to a partially        dependent Dirichlet process pdDP defined as

$\begin{matrix}{\mspace{79mu} {{{H( \cdot )} = {\sum\limits_{k = 1}^{\infty}\; {w_{k}\text{?}\text{?}( \cdot )}}}{\text{?}\text{indicates text missing or illegible when filed}}}} & (3)\end{matrix}$

-   -   δ_(ξ)(·) the Dirac measurement at ξ.

N(m, Σ) designates multi-variate normal-Gaussian law in dimension p=3parameterized by its mean m and its covariance matrix Σ, of density:

$\mspace{79mu} {{f_{}(x)} = {\frac{1}{\left( {2\; \pi} \right)\text{?}{\sum }^{\frac{1}{2}}}{\exp\left( {{- \frac{1}{2}}\left( {x - m} \right)^{T}{\sum\limits_{\;}^{- 1}\; \left( {x - m} \right)}} \right)}}}$?indicates text missing or illegible when filed

The generative model defined in the expression (2) and the concept ofpdDP of the expression (3) are explained hereinbelow.

The expression of H in (3) is a statistical mixture in which w_(k)corresponds to the weight of the K^(th) component, {tilde over (Z)}_(k)to the spatial parameters and {tilde over (Q)}_(k) to the kinetic, thatis to say the temporal probability distribution, of the K^(th)space-time component of the mixture. The distributions for the w_(k) and{tilde over (Z)}_(k) are identical to those of a Dirichlet process.

w=w₁, w₂, . . . ˜GEM(θ) applies, as indicated in the paper by J. Pitmanentitled “Some developments of the Blackwell-MacQueen urn scheme”, T. F.et al., editor, Statistics Probability and Game Theory, Papers in honorof David Blackwell, volume 30 of Lecture Notes-Monograph Series, pages245-267, Institute of Mathematical Statistics, 1996, such that:

-   -   w₁=v₁;    -   for any k≧2, w_(k)=v_(k)·Π_(i−1) ^(k−1)(1 v_(i)) in which v₁,        v₂, . . . are identically distributed according to a law Beta        (1, θ), of which the can be obtained by using the following        expression:

${f_{Beta}\left( {{\upsilon a},b} \right)} = {\frac{\Gamma \left( {a + b} \right)}{{\Gamma (a)}{\Gamma (a)}}{\upsilon^{a - 1}\left( {1 - \upsilon} \right)}^{b - 1}}$

in which the symbol Γ represents the gamma function.

{tilde over (Z)}={tilde over (Z)}₁, {tilde over (Z)}₂, . . . areidentically distributed according to the a priori law G₀ which can be,for example, chosen as an inverse normal-Wishart model (NIW_(ρ, n) ₀_(,μ) ₀ _(,Σ) ₀ ) such that {tilde over (Z)}_(k)=(m_(k), Σ_(k)) with:

m _(k)|Σ_(k)

N(m _(k)|μ₀, Σ_(k)/ρ)

Σ_(k) ¹

W(Σ_(k) ¹ |n ₀,(n ₀Σ₀)¹).

in which W designates the Wishart distribution for which the density isexpressed by using the expression:

$\begin{matrix}{\mspace{79mu} {{{f_{}\left( \sum\limits_{\;}^{- 1} \right)} = {\frac{{{n_{0}\sum\limits_{0\;}}}\text{?}}{\text{?}{\sum\limits_{\;}}\text{?}}{\exp\left( {{- \frac{n_{0}}{2}}{{Tr}\left( {\underset{\;}{\sum\limits_{\;}^{- 1}}\sum\limits_{0\;}} \right)}} \right)}}}{\text{?}\text{indicates text missing or illegible when filed}}}} & (4)\end{matrix}$

in which n₀ is the degree of freedom of the Wishart law and Σ₀ theassociated scale matrix.

This is where the spatial concentration becomes involved as endogenousco-variable for the distribution of the kinetic parameters {tilde over(Q)}_(k). First of all, the spatial distribution S(·) is introduced asmarginal law of the space-time distribution H. The expression of S(·) isgiven by:

$\mspace{79mu} {{S( \cdot )} = {{\int_{\mathcal{M}{()}}^{\;}{H\left( {\cdot {,{}}} \right)}} = {\sum\limits_{k = 1}^{\infty}\; {w_{k}\text{?}( \cdot )}}}}$?indicates text missing or illegible when filed

in which M(T) represents the set of the measurements over T.

It can be noted that, by the definitions of w_(k) and {tilde over(Z)}_(k), S(·) is distributed according to a Dirichlet process ofconcentration parameter 0 and a basis measurement G₀:

S˜DP(θ, G₀)

{tilde over (Q)}={tilde over (Q)}₁, {tilde over (Q)}₂. . . are thenindependently distributed according to:

{tilde over (Q)} _(k) |{tilde over (Z)} _(k) , S

K ₀ ^(λ) ^(k)

with λ_(k)=φ_(λ)({tilde over (Z)}_(k), S) in which φ_(λ)(·) takes itsvalues from Λ⊂

. For example, the function φ_(λ) can be chosen as being the value ofthe expectation of the spatial density taken over the component K:

  λ_(k) = ?(?w_(m)(??))   with      ? ∼ (??)?indicates text missing or illegible when filed

It can be noted that, in the case where {tilde over (Q)}_(k) isindependent of the marginal S such that λ_(k)=φ*_(λ)({tilde over(Z)}_(k)), we obtain:

{tilde over (Q)} _(k) |{tilde over (Z)} _(k)

K ₀ ^(φ*) ^(k) ^(({tilde over (Z)}) ^(k) ⁾

{tilde over (Z)} _(k) , {tilde over (Q)} _(k)

F ₀

expressions in which F₀({tilde over (Z)}_(k), {tilde over(Q)}_(k))=G₀({tilde over (Z)}_(k))K₀ ^(φ*) ^(k) ^(({tilde over (Z)})^(k) ⁾({tilde over (Q)}_(k)). This implies that H˜DP(θ, F₀).

The distribution of H is a Dirichlet process. In this situation theconcentration is not involved as endogenous co-variable. However, themodel then covers the possibility of having recourse to an exogenousco-variable originating from information of anatomical type for example.

The distribution of the K₀ ^(λ) remains to be established. K₀ ^(Λ)={K₀^(λ): λ∈Λ} is defined as an innumerable collection of random probabilitymeasurements which accepts, as distribution, a local Dirichlet processdenoted IDP, defined as follows. If:

-   V=V₁, V₂, . . . are identically distributed such that V_(j)    Beta (1, v),-   Q*=Q*₁, Q*₂, . . . are identically distributed such that Q*_(j)    PT_(L(A, Q) ₀) is a Pólya tree with L-levels, of parameters A and of    basis distribution Q₀,-   Γ=Γ₁, Γ₂, . . . are identically distributed such that Γ_(j)    γ₀,    then,

$\mspace{20mu} {{K_{0}^{\lambda}( \cdot )} = {\sum\limits_{j = 1}^{\infty}\; {{p_{j}(\lambda)}\text{?}( \cdot )}}}$?indicates text missing or illegible when filed

with p(λ)=p₁(λ),p₂(λ), . . . defined by the expression in which:

${p_{j}(\lambda)} = {{{\overset{\sim}{V}}_{j}(\lambda)}{\prod\limits_{l < j}\; \left( {1 - {{\overset{\sim}{V}}_{l}(\lambda)}} \right)}}$in  which${{\overset{\sim}{V}}_{j}(\lambda)} = {V_{j}\left( {{d\left( {\lambda,\Gamma_{j}} \right)} < \psi} \right)}$

and d(λ, λ′) represents the Euclidian distance between λ and λ′.

It will be noted that if K₀ ^(Λ)−{K₀ ^(λ): λ∈Λ}˜IDP(v, PT(A, Q₀), γ₀,ψ), then for any λ∈Λ, K₀ ^(λ)|λ˜DP(v, PT(A, Q₀)). The conditional law ofK₀ ^(λ) knows that λ is a Pólya tree Dirichlet process.

It will be noted that in the limit case where ψ→∞, there is, for anyλ∈Λ, K₀ ^(λ)=K₀ with K₀˜DP(v, PT(A, Q₀)), consequently

H˜DP(θ, G ₀ ×K ₀)

The definition of a Pólyatree to complete the model remains to bereviewed. We adopt the Lavine formulism presented in his paper entitledSome aspects of Pólya tree distributions for statistical modelling, Ann.Statist., 20:1222-1235, 1992, for the definition of the dyadic Pólyatree. Take E={0, 1} and the E^(l) the l-cartesian product E× . . . ×Ewith E⁰=. Let E*=∪_(l=0) ^(∞)E^(l). Let B_()={T} and let Π={B₀, B₁,B₀₀, B₀₁, . . . } and separable binary tree of partitions of T such thatsets B_(∈) ₁ _(. . . ∈l+1) are obtained by binary subdivision of thesets B_(∈) ₁ _(. . . ∈l) and ∪_(l=0) ^(∞)B_(∈) ₁ _(. . . ∈l+1) definesthe measurable sets. The index 0 (respectively 1) is interpreted as 0(respectively 1).

It is said that the random measurement Q over T accepts as distributiona Pólya tree of parameters (A, Q₀) that will be noted Q˜PT(A, Q₀), ifthere are positive real values A={α_(∈): ∈∈E*} and random variablesW={W_(∈): ∈∈E*} such that:

-   -   the collection W consists of mutually independent random        variables,    -   ∀_(∈)∈E*, W_(∈)˜Beta (α_(∈0), α_(∈1)),    -   ∀l=1, 2, . . . et ∀_(∈)∈E^(l),

${Q\left( B_{\varepsilon_{1}\; \ldots \mspace{11mu} \varepsilon_{l}} \right)} = {\left( {\prod\limits_{\underset{\{{{k:\varepsilon_{k}} = 0}\}}{k = 1}}^{l}\; W_{\varepsilon_{1}\mspace{11mu} \ldots \mspace{11mu} \varepsilon_{k - 1}}} \right)\left( {\prod\limits_{\underset{\{{{k:\varepsilon_{k}} = 1}\}}{k = 1}}^{l}\; W_{\varepsilon_{1}\mspace{11mu} \ldots \mspace{11mu} \varepsilon_{k - 1}}} \right)}$

with factors equal to W_() in which 1−W_() if k−1.Consider the probability measurement Q₀. The Pólya tree can be centeredaround Q₀ by taking:

$B_{\varepsilon_{1}\; \ldots \mspace{11mu} \varepsilon_{l}} = \left( {{Q_{0}^{- 1}\left( \frac{\sum\limits_{k = 1}^{l}{\varepsilon_{k}2^{k - 1}}}{2^{l}} \right)},{Q_{0}^{- 1}\left( \frac{{\sum\limits_{k = 1}^{l}{\varepsilon_{k}2^{k - 1}}} + 1}{2^{l}} \right)}} \right)$and  α_(ε 0) = α_(ε 1), ∀ε ∈ E^(*)

In this case, we note α_(∈)=a_(l) for l=1, 2, . . . and ∀_(∈)∈E^(l).This construction scheme defines the canonical Pólya tree. In thiscontext, if Q˜PT(A, Q₀), then

(Q)=Q₀.

The construction of the partition can be continued to infinity but,since it is not possible to consider computing an infinite tree inpractice, the Pólya trees are often simplified by specifying theparameters and by constructing the partition up to a bounded level L. Ifwe note f_(Q)(·), the density of probability of Q˜PT_(L)(A, Q₀), then:

∀t∈B _(∈) ₁ _(. . . ∈) _(L) , f _(Q)(t)=2^(L) f _(Q) _(n) (t)·Q(B _(∈) ₁_(. . . ∈) _(L) )

The justification for a partial specification of the Pólya tree can besupported by the fact that, for the continuous distributions, theparameters rapidly increase as a function of the levels of the tree. Theupdating of the law from the data then has little a priori influence onthe law because the number of observations in the high levels becomesmall compared to the a priori parameter.

Since the objective is to culminate with an algorithm for reconstructionof the functional volumes, the method according to the invention can beprovided with auxiliary variables which will render the iterativereconstruction procedure and the implementation of an MCMC-type method,MCMC standing for “Markov chain Monte-Carlo”, effective.

We introduce C=C₁, C₂, . . . , C_(n), the data classification variables,to the components of H such that (Z_(i), Q_(i))=({tilde over (Z)}_(C)_(i) , {tilde over (Q)}_(C) _(i) ) for any i<n. Similarly, let D=D₁, D₂,. . . such that {tilde over (Q)}_(k)=Q*_(D) _(k) for any K.

The auxiliary variables for the sampling by slices, usually referred toby the expression “slice sampling”, are described hereinbelow.

In order to implement an effective Gibbs sampler, we follow a so-calledslice sampling technique as described in the paper by M. Kalli, J. E.Griffin and S. G. Walker entitled Slice sampling mixture models,Statistics and Computing, 21(1):93-105, 2011. This approach requires theintroduction of additional variables.

Let u=u₁, u₂, . . . , u_(n) be uniform random variables. The jointdensity of (X_(i), T_(i), u_(i)), for any positive sequence (ξ_(k)), canbe determined by using the expression:

$\left( {X_{i},T_{i},\left. u_{i} \middle| w \right.,\overset{\sim}{Z},\overset{\sim}{Q}} \right) \sim {\sum\limits_{k = 1}^{\infty}{{\left( {\left. u_{i} \middle| 0 \right.,\xi_{k}} \right)}w_{k}{\left( X_{i} \middle| {\overset{\sim}{Z}}_{k} \right)}{{\overset{\sim}{Q}}_{k}\left( T_{i} \right)}}}$

-   in which u(·|a, b) is the uniform distribution over ]a, b]. We    easily check that, by marginalization on u_(i) of the preceding    formula, we obtain the distribution of (X_(i), T_(i)|w, {tilde over    (Z)}, {tilde over (Q)}). We propose adopting a dependent slice    variable (ξ_(k)), that is to say that, for any K, ξ_(k)=min (w_(k),    ζ)-   in which ζ∈]0, 1] and is independent of w_(k). In practice, values    of ζ corresponding to the expectation of the weight of the first    component that does not have any datum allocated offer satisfactory    mixture properties

$\left( {\zeta \approx \frac{0}{\left( {n + 0} \right)\left( {1 + 0} \right)}} \right).$

-    Provided with this choice of ξ_(k), the following expression is    obtained:

$\left( {X_{i},T_{i},\left. u_{i} \middle| w \right.,\overset{\sim}{Z},\overset{\sim}{Q}} \right) \sim {{\zeta^{- 1}{\sum\limits_{u_{i} < \zeta < w_{k}}{w_{k}{\left( X_{i} \middle| {\overset{\sim}{Z}}_{k} \right)}{{\overset{\sim}{Q}}_{k}\left( T_{i} \right)}}}} + {\sum\limits_{u_{i} < w_{k} \leq \zeta}{{\left( X_{i} \middle| {\overset{\sim}{Z}}_{k} \right)}{{\overset{\sim}{Q}}_{k}\left( T_{i} \right)}}}}$

-   in which the two sums are finite by virtue of #{j: w_(j)>ε}<∞ for    any ε>0.

Similarly, let v=v₁, v₂, . . . be uniform auxiliary variables and ζ∈]0,1]. Then, for any K, the following applies:

$\left( {{\overset{\sim}{Q}}_{k},\left. \upsilon_{k} \middle| V \right.,Q^{*},\Gamma,w,{\overset{\sim}{Z}}_{k},S} \right) \sim {{\varsigma^{- 1}{\sum\limits_{\upsilon_{k} < \varsigma < {p_{j}{(\lambda_{k})}}}{{p_{j}\left( \lambda_{k} \right)}{\delta_{Q_{j}^{*}}\left( {\overset{\sim}{Q}}_{k} \right)}}}} + {\sum\limits_{\upsilon_{k} < {p_{j}{(\lambda_{k})}} \leq \varsigma}{\delta_{Q_{j}^{*}}\left( {\overset{\sim}{Q}}_{k} \right)}}}$

expression in which λ_(k)=φ_(λ)({tilde over (Z)}_(k), S) andp_(j)(λ_(k)) are defined above.

The model (2) is then rewritten by adding to it the hidden and auxiliaryvariables:

$\begin{matrix}{\mspace{79mu} {{\left. Y_{i} \middle| {X_{i}\overset{ind}{\sim}{\left( Y_{i} \middle| X_{i} \right)}} \right.\mspace{20mu} {X_{i},\left. T_{i} \middle| C_{i} \right.,{D\overset{ind}{\sim}{{\left( X_{i} \middle| {\overset{\sim}{Z}}_{C_{i}} \right)} \times {Q_{D_{C_{i}}}^{*}\left( T_{i} \right)}}}}}\mspace{20mu} {C_{i},\left. u_{i} \middle| {w\overset{iid}{\sim}{\sum\limits_{k = 1}^{\infty}{{\left( {\left. u_{i} \middle| 0 \right.,{\min \left( {\zeta,w_{k}} \right)}} \right)}w_{k}{\delta_{k}\left( C_{i} \right)}}}} \right.}{D_{k},\left. \upsilon_{k} \middle| V \right.,\Gamma,w,{\overset{\sim}{Z}\overset{ind}{\sim}{\sum\limits_{j = 1}^{\infty}{{\left( {\left. \upsilon_{k} \middle| 0 \right.,{\min \left( {\varsigma,{p_{j}\left( \lambda_{k} \right)}} \right)}} \right)}{p_{j}\left( \lambda_{k} \right)}{\delta_{j}\left( D_{k} \right)}}}}}\mspace{20mu} {w \sim {{GEM}(\theta)}}\mspace{20mu} {{\overset{\sim}{Z}}_{k}\overset{iid}{\sim}{\mathcal{I}}_{\varrho,n_{0},\mu_{0},\Sigma_{0}}}\mspace{20mu} {V_{j}\overset{iid}{\sim}{{Beta}\left( {1,\nu} \right)}}\mspace{20mu} {\Gamma_{j}\overset{iid}{\sim}\mathrm{\Upsilon}_{0}}\mspace{20mu} {Q_{j}^{*}\overset{iid}{\sim}{{PT}_{L}\left( {,Q_{0}} \right)}}\mspace{20mu} {with}\mspace{20mu} {{{p_{j}\left( \lambda_{k} \right)} = {{{\overset{\sim}{V}}_{j}\left( \lambda_{k} \right)}{\Pi_{l < j}\left( {1 - {{\overset{\sim}{V}}_{l}\left( \lambda_{k} \right)}} \right)}}},\mspace{20mu} {{{\overset{\sim}{V}}_{j}\left( \lambda_{k} \right)} - {V_{j}\left( {{d\left( {\lambda_{k},\Gamma_{j}} \right)} < \psi} \right)}}}\mspace{79mu} {and}\mspace{79mu} {\lambda_{k} = {_{{\overset{\sim}{x}}_{k}}\left( {\sum\limits_{l = 1}^{\infty}{w_{l}{\left( {\overset{\sim}{x}}_{k} \middle| {\overset{\sim}{Z}}_{l} \right)}}} \right)}}}} & (5)\end{matrix}$

From this generative model, it is possible to express the different aposteriori conditional laws which will be involved in the Gibbs sampler.

The sampling of the a posteriori conditional laws by Markov chainMonte-Carlo MCMC algorithm is described hereinbelow with the support ofFIGS. 2 to 7.

From the expression (5), the MCMC sampler generates, in succession,samples originating from the following conditional distributions:

-   Places of emission/allocation to the components of H: (X,    C|w,u,{hacek over (Z)}, Q*, D, Y, T)-   Spatial parameters: ({hacek over (Z)}|C, X)-   Weights and auxiliary variables of H: (w, u|C)-   Allocation of the kinetics to the components of K₀ ^(Λ): (D|V, v,    Q*, Γ, w, {hacek over (Z)}, C, T)-   Pólya trees of the kinetics: (Q*|D, C, T)

Weightings, local variables and auxiliaries of K₀ ^(Λ): (V, Γ, v|D, w,{hacek over (Z)}, T)

FIG. 2 illustrates how the weights and auxiliary variables of thedistribution H of the space-time components are generated, that is tosay how samples of conditional distribution (w,u|C) are generated.

We denote k_(n) the latent number of components of the space-timemixture for n observations and n_(k)=#{i: C_(i)=k} for k≦k_(n). We adoptthe approach presented in the paper by Kalli et al., previously cited,and proceed with the joint generation of (w, u|C) by first of allsampling w₁, w₂, . . . , w_(k) _(n) |C, then u|w₁, w₂, . . . , w_(k)_(n) , C, and finally w_(k) _(n) ₊₁, w_(k) _(n) ₊₂, . . . |u, w₁, w₂, .. . , w_(k) _(n), C.

In this context, by following the new definition and the corollarytwenty of the paper by J. Pitman entitled Some developments of theBlackwell-MacQueen urn scheme in T. F. et al., editor, Statistics,Probability and Game Theory; Papers in honor of David Blackwell, volume30 of Lecture Notes—Monograph Series, pages 245-267, Institute ofMathematical Statistics, 1996, we are able to express the a posteriorilaw of the Dirichlet process as the sum of a finite sum of componentswhose weights are distributed by a Dirichlet distribution and of adata-independent Dirichlet process. The presence of the auxiliaryvariable u makes it possible to generate only a finite number k* ofweights and of places of emission in the generation of the independentDirichlet process.

In a first stage 200, the w_(k) are generated for k≦k_(n).

(w ₁ , w ₂ , . . . , w _(k) _(n) , r _(k) _(n) |C)˜Dirichlet (n ₁ , n ₂, . . . , n _(k) _(n) , θ)

In a second stage 201, the samples u_(i)|w₁, w₂, . . . , w_(k) _(n) , Care generated.

u _(i) |C

u(u _(i)|0, min (w _(C) _(i) , ζ))

u*=min(u)

In a third stage 202, the samples w_(k) are generated for k>k_(n).As long as r_(k−1)>u*,

U _(k)˜Beta (1, θ)

w _(k) =U _(k) r _(k−1)

r _(k) =r _(k−1)(1−U _(k))

In a fourth stage 203, k*=min({k: r_(k)<u*}) is posited.Clearly, w_(k)<u* for any k>k*. For this reason, a finite set of w_(k)is sufficient for the allocation of the data to the components. Itshould be added that we also abandon in the computations all the w_(k)such that w_(k)<u* for k≦k_(n).

FIG. 3 schematically illustrates how the spatial parameters aregenerated, that is to say how samples of conditional distribution {tildeover (Z)} C, X) are generated.

In a first stage 300, ({tilde over (Z)}_(k)|C, X) are generated for anyk≦k_(n), from an inverse normal-Wishart law whose parameters are updatedfrom the observations, through the conjugation property:

$\left. m_{k} \middle| \Sigma_{k} \right.,C,{X\overset{ind}{\sim}{\left( {\mu_{k},\frac{\Sigma_{k}}{\varrho_{k}}} \right)}}$with 𝜚_(k) = 𝜚 + n_(k)$\mu_{k} = \frac{{\varrho \; \mu_{0}} + {n_{k}}}{\varrho + n_{k}}$$= {\frac{1}{n_{k}}{\sum\limits_{\{{{i:C_{i}} = k}\}}X_{i}}}$ and$\left. \Sigma_{k}^{- 1} \middle| C \right.,{X\overset{ind}{\sim}{\left( {\eta_{k},\left( {\eta_{k}{\overset{\sim}{\Sigma}}_{k}} \right)^{- 1}} \right)}}$with η_(k) = n₀ − n_(k)${\overset{\sim}{\Sigma}}_{k} = {\frac{1}{\eta_{k}}\left\lbrack {{n_{0}\Sigma_{0}} + {n_{k}} + \frac{\left( {\mu_{0} -} \right)\left( {\mu_{0} -} \right)^{\top}}{\frac{1}{\varrho} + \frac{1}{n_{k}}}} \right\rbrack}$$= {\frac{1}{n_{k}}{\sum\limits_{\{{{i:C_{i}} = k}\}}{\left( {X_{i} -} \right)\left( {X_{i} -} \right)^{\top}}}}$

In a second stage 301, the samples {tilde over (Z)}_(k) are generated:

{tilde over (Z)}_(k)

NIW_(ρ, n) ₀ _(,μ) ₀ _(,Σ) ₀ , for k_(n)<k≦k*

FIG. 4 schematically illustrates one way of determining the local andauxiliary variables of the set K₀ ^(λ).

For this, (V, Γ, v|D, w, {tilde over (Z)}, T) must be generated.

We consider as a priori law γ₀(·)=u(·|a_(Γ), b_(Γ)). For any j, wedefine

D _(j) ={k≦k _(n) : D _(k) =j}

D _(j) ⁺ ={k≦k _(n) : D _(k) >j}

and designate, for any j, m_(j)=#D_(j) as the number of components of Hallocated to the j^(th) component of K₀ ^(Λ).We also define the greatest index of the components of K₀ ^(Λ) having atleast one component of H allocated:

J _(k) _(n) =max({j: m _(j)>0})

We adopt the same generation scheme as for the sampling of (w, u|C),that is to say we first generate (V_(j), Γ_(j)|D, w, {tilde over (Z)},T) for j≦J_(k) _(n) 400, then

$\left( {{vD},V_{1},\cdots \mspace{14mu},V_{J_{\kappa_{n}}},\Gamma_{1},\cdots \mspace{14mu},\Gamma_{J_{\kappa_{n}}},w,\overset{\sim}{Z}} \right)$

401 and, finally, V_(j) and Γ_(j) for j>j_(k) _(n) 402.For any j≦J_(k) _(n) , the following are generated independently:

$\left. V_{j} \middle| D \right.,{T\overset{ind}{\sim}{{Beta}\left( {{1 + m_{j}},{\nu {\sum\limits_{l = {j + 1}}^{J_{\kappa_{n}}}m_{l}}}} \right)}}$

For any j≦J_(k) _(n) , we posit:

$a_{j,\Gamma}^{\prime} = {\max\left( {{\max\limits_{k \in _{j}}\left( {\lambda_{k} - \psi} \right)},a_{\Gamma}} \right)}$$b_{j,\Gamma}^{\prime} = {\min\left( {{\min\limits_{k \in _{j}}\left( {\lambda_{k} + \psi} \right)},b_{\Gamma}} \right)}$

then generate:

$\left. \Gamma_{j} \middle| V_{j} \right.,D,w,{\overset{\sim}{Z}\overset{ind}{\sim}{_{j}^{*}\left( \Gamma_{j} \right)}}$in  which${_{j}^{*}\left( \Gamma_{j} \right)} \propto {\left( {a_{j,\Gamma}^{\prime} < \Gamma_{j} < b_{j,\Gamma}^{\prime}} \right) \times {\prod\limits_{k \in _{j}^{+}}\; \left( {1 - {V_{j}\left( {{d\left( {\lambda_{k},\Gamma_{j}} \right)} < \psi} \right)}} \right)}}$

The sampling of u*_(j) requires a Metropolis-Hastings (MH) step. Wepropose an independent MH sampler with, as candidate distribution:

g(Γ_(j))=u(Γ_(j) |a′ _(j,Γ) , b′ _(j,Γ))

Thus, on the iteration \|⊥, for any j≦J_(k) _(n),

-   -   generate Γ*_(j)        u(Γ_(j)|a′_(j,Γ), b′_(j,Γ)).    -   Compute the acceptance ratio w(Γ*_(j), Γ_(j) ^((t)))

${\omega \left( {\Gamma_{j}^{*} \cdot \Gamma_{j}^{(t)}} \right)} = \begin{matrix}{\prod\limits_{k \in _{j}^{+}}\; \left( {1 - {V_{j}\left( {{d\left( {\lambda_{k},\Gamma_{j}^{*}} \right)} < \psi} \right)}} \right)} \\{\prod\limits_{k \in _{j}^{+}}\; \left( {1 - {V_{j}\left( {{d\left( {\lambda_{k},\Gamma_{j}^{(t)}} \right)} < \psi} \right)}} \right)}\end{matrix}$

in which Γ_(j) ^((t)) is the value of the Markov chain on the iterationt.

-   -   Test the acceptance of the proposition Γ*_(j), i.e. assign Γ_(j)        ^((t+1))=Γ*_(j) with a probability

min(1, ω(Γ*_(j), Γ_(j) ^((t)))),

otherwise assign Γ_(j) ^((t+1))=Γ_(j) ^((t)).

-   -   For any k≦k_(n), posit ρ_(j,k)=ρ_(j−1,k)(1−V_(j)        (d(λ_(k), Γ_(j))<ψ)) and ρ*_(j)=min(ρ_(j,1), . . . , ρ_(j,k)        _(n) ).    -   For any k≦k_(n), independently generate

$\left. \upsilon_{k} \middle| D \right.,V_{1},\ldots \mspace{14mu},V_{J_{\kappa_{n}}},\Gamma_{1},\ldots \mspace{14mu},\Gamma_{J_{\kappa_{n}}},w,{\overset{\sim}{Z}\overset{ind}{\sim}{\left( {\left. \upsilon_{k} \middle| 0 \right.,{\min \left( {{p_{D_{k}}\left( \lambda_{k} \right)},\varsigma} \right)}} \right)}}$

and posit

v*=min(v)

-   -   Generate V_(j) et Γ_(j) for j>J_(k) _(n) . As long as        ρ*_(j−1)>v*,

V _(j)

Beta(1, v)

Γ_(j)

u(Γ_(j) |a _(Γ) , b _(Γ))

-   -   For any k≦k_(n), posit ρ_(j,k)=ρ_(j−1,k)(1−V_(j)        (d(λ_(k), Γ_(j))<ψ)) and ρ*_(j)=min(ρ_(j,1), . . . , ρ_(j,k)        _(n) ).    -   Posit J*=min({j: ρ*_(j)<v*}).

FIG. 5 illustrates how the kinetics are allocated to the components ofK₀ ^(λ), that is to say how (D|V, v, Q*, Γ, w, {tilde over (Z)}, T) isgenerated.

In a first stage 500, the following is generated independently for anyk≦k_(n):

$\left( {\left. D_{k} \middle| V \right.,v,Q^{*},\Gamma,w,\overset{\sim}{Z},C,T} \right)\overset{ind}{\sim}{\sum\limits_{j = 1}^{J^{*}}{p_{j,k}{\delta_{j}( \cdot )}}}$with${p_{j,k} \propto {\left( {{p_{j}\left( \lambda_{k} \right)} > \upsilon_{k}} \right){\max \left( {{p_{j}\left( \lambda_{\kappa} \right)},\varsigma} \right)}{\prod\limits_{\{{{i:C_{i}} = k}\}}\; {f_{Q_{j}^{*}}\left( T_{i} \right)}}}}$

and Σ_(j=1) ^(J*)P_(j,k)−1.In a second stage 501, the following is generated for any K, such thatk_(n)<k<k*:

$\left. D_{k} \middle| V \right.,v,\Gamma,w,{\overset{\sim}{Z}\overset{ind}{\sim}{\sum\limits_{j = 1}^{J^{*}}{p_{j,k}{\delta_{j}( \cdot )}}}}$with p_(j, k) ∝ (p_(j)(λ_(k)) > υ_(k))max (p_(j)(λ_(κ)), 𝜍)${{and}\mspace{14mu} {\sum\limits_{j = 1}^{J^{*}}p_{j,k}}} - 1.$

The Pólya trees of the generated kinetics must then be generated.According to the conjugation property of the Pólya tree, for any j≦J_(k)_(n) , (Q*_(j)|D, C, T) is generated independently according to a Pólyatree distribution with updated parameters A′_(j)−{α′_(j,τ): ε∈E*} suchthat

Q* _(j) |D, C, T

PT _(L)(A′ _(j) , Q ₀)

with

α′_(j,τ)−α_(τ)+n_(j,τ)

and

n _(j,∈) =#{i∈{1, . . . , n}: (D _(C) _(i) =j)̂(T _(i) ∈B _(∈))}

Then, for any J_(k) _(n) <j≦J*, independently generate

Q* _(j) |D, C, T

PT _(L)(A, Q ₀)

FIG. 6 shows how the places of emission are determined and how thecomponents of H are determined. This is by generating (C, X|w, u, {tildeover (Z)}, Q*, D, Y, T).

The conditional law has a form which does not allow for a directgeneration. In effect, it involves a projector P(Y_(i)|X_(i)) whichincorporates geometrical and physical characteristics of the detectiondevice, which in practice precludes finding an analytical expression forthe distribution that offers a hope of simple generation. We thereforehave recourse to a Metropolis-Hastings step. We begin by expressing theconditional distribution by using the Bayes rule:

$C_{i},\left. X_{i} \middle| Y_{i} \right.,T_{i},u_{i},w,\overset{\sim}{Z},Q^{*},{D\overset{ind}{\sim}{{\left( Y_{i} \middle| X_{i} \right)} \times {{\overset{\sim}{G}}_{i}\left( {C_{i},X_{i}} \right)}}}$with${f_{{\overset{\sim}{G}}_{i}}\left( {C_{i},X_{i}} \right)} \propto {\sum\limits_{k = 1}^{\kappa^{*}}{\left( {w_{k} > u_{i}} \right){\max \left( {w_{k},\zeta} \right)}{f_{}\left( X_{i} \middle| {\overset{\sim}{Z}}_{k} \right)}{f_{Q_{D_{k}}^{*}}\left( T_{i} \right)}{\delta_{k}\left( C_{i} \right)}}}$

and P(Y_(i)|X_(i)) corresponds to the likelihood of the projector forthe sets of observed events (truncated). Since this likelihood does notaccept any standard form allowing for direct sampling, we use anindependent Metropolis-Hastings algorithm and have to define a candidatelaw that we denote G*_(i) for i=1, . . . , n. On each iteration t+1 wehave to generate (C*_(i), X*_(i))˜G*_(i) and compute the acceptanceratio. This step is the most costly in terms of computations.

We thus propose a candidate law of the form:

f _(G*) _(i) (C* _(i) , X* _(i))∝f _(N)(X* _(i) |m* _(Y) _(i) , Σ*_(Y)_(i) )×f _(G) _(i) (C* _(i) , X* _(i))

We replace the projector P(Y_(i)|X_(i)) with N(X_(i)|m*_(Y) _(i) ,Σ*_(Y) _(i) ) in the conditional law to construct the candidatedistribution 600, in which m*_(Y) _(i) and Σ*_(Y) _(i) are chosen insuch a way as to suitably approximate P(Y_(i)|X_(i)). We can thenexpress f*_(G) _(i) (·) by using the expression:

${f_{G_{i}^{*}}\left( {C_{i}^{*},X_{i}^{*}} \right)} \propto {\sum\limits_{k = 1}^{\kappa^{*}}{\left( {w_{k} > u_{i}} \right){\max \left( {w_{k},\zeta} \right)}{\overset{\_}{\eta}}_{k,Y_{i}}{f_{}\left( {\left. X_{i}^{*} \middle| {\overset{\_}{m}}_{k,Y_{i}} \right.,{\overset{\_}{\Sigma}}_{k,Y_{i}}} \right)}{f_{Q_{D_{k}}^{*}}\left( T_{i} \right)}{\delta_{k}\left( C_{i}^{*} \right)}}}$

with, for any k≦k* and m≦M,

η _(k,m) =f _(N)(m* _(m) |m _(k), Σ*_(m)−Σ_(k))

m _(k,m)=((Σ*_(m))⁻¹+Σ_(k) ⁻¹)⁻¹((Σ*_(m))⁻¹ m* _(m)+Σ_(k) ⁻¹ m _(k))

Σ _(k,m)=((Σ*_(m))⁻¹+Σ_(k) ⁻¹)⁻¹

It is then possible to proceed with the sampling of C_(i), X_(i)), forany i≦n on the iteration t+1.C*_(i) is generated independently 601 such that:

$C_{i}^{*}\overset{ind}{\sim}{\sum\limits_{k = 1}^{\kappa^{*}}{{\overset{\sim}{w}}_{k,i}{\delta_{k}( \cdot )}}}$with${\overset{\sim}{w}}_{k,i} \propto {\left( {w_{k} > u_{i}} \right){\max \left( {w_{k},\zeta} \right)}{\overset{\sim}{\eta}}_{k,Y,}{f_{Q_{D_{k}}^{*}}\left( T_{i} \right)}}$and ${\sum\limits_{k = 1}^{\kappa^{*}}{\overset{\_}{w}}_{k,i}} = 1.$

Then X*_(i) is generated independently 602 such that:

X* _(i)

N(X* _(i) | m _(C*) _(i) _(,Y) _(i) , Σ _(C*) _(i) _(,Y) _(i) )

The acceptance ratios ω((C*_(i), X*_(i)), (C_(i) ^((t)), X_(i) ^((t))))are then computed 603:

$\begin{matrix}{{\omega \left( {\left( {C_{i}^{*},X_{i}^{*}} \right),\left( {C_{i}^{(t)},X_{i}^{(t)}} \right)} \right)} = {\frac{\left( Y_{i} \middle| X_{i}^{*} \right)}{\left( Y_{i} \middle| X_{i}^{(t)} \right)} \cdot \frac{f_{}\left( {\left. X_{i}^{(t)} \middle| m_{Y_{i}}^{*} \right.,\Sigma_{Y_{i}}^{*}} \right)}{f_{}\left( {\left. X_{i}^{*} \middle| m_{Y_{i}}^{*} \right.,\Sigma_{Y_{i}}^{*}} \right)}}} \\{= {\omega \left( {X_{i}^{*},X_{i}^{(t)}} \right)}}\end{matrix}$

It will be able to be noted that {tilde over (G)}_(i)(C_(i), X_(i)) andC_(i) disappear from the acceptance ratio.

It then remains to test 604 the acceptance of (C*_(i), X*_(i)), that isto say assign (C_(i) ^((t+1)), X_(i) ^((t+1)))=(C*_(i), X*_(i)) with theprobability

min(1,ω(X*_(i), X_(i) ^((t))))

otherwise assign (C_(i) ^((t+1)), X_(i) ^((t+1)))=(C_(i) ^((t)), X_(i)^((t))))FIG. 7 presents a way of computing the concentration λ_(k). Theevaluation of λ_(k) is necessary for the generation of (V, Γ, v|D, w,{tilde over (Z)}, T) and of (D|V, v, Q*, Γ, w, {tilde over (Z)}, C, T).The expression of λ_(k) for the particular choice of φ_(λ)({tilde over(Z)}_(k), S) given previously, can be determined according to theequation (5) by using the expression:

$\lambda_{k} = {_{{\overset{\sim}{x}}_{k}}\left( {\sum\limits_{m = 1}^{\infty}{w_{m}{\left( {\overset{\sim}{x}}_{k} \middle| {\overset{\sim}{Z}}_{m} \right)}}} \right)}$

This expression can be evaluated by a significant sampling technique.For this, for a given number of particles T_(λ) and for any k, thefollowing is generated independently 700 for any l≦T_(λ):

{tilde over (x)} _(kl)

N({tilde over (x)} _(kl) |{tilde over (Z)} _(k))

Then, the following is computed 701 for any l≦T_(λ):

${\overset{\sim}{\lambda}}_{kl} \approx {\sum\limits_{m = 1}^{\kappa^{*}}{w_{m}{\left( {\overset{\sim}{x}}_{kl} \middle| {\overset{\sim}{Z}}_{m} \right)}}}$

In a third stage 702, λ_(k) is evaluated:

$\lambda_{k} \approx {\frac{1}{T_{\lambda}}{\sum\limits_{l = 1}^{T_{\lambda}}{\overset{\sim}{\lambda}}_{kl}}}$

The functional volumes VF then remain to be determined 105.

Once the MCMC algorithm has reached balance equilibrium, the samplerprovides draws of the a posteriori joint law of (X, C, D, w, u, {tildeover (Z)}, Q*|Y, T). For each draw of parameters retained, it ispossible to construct a sample of the space-time activity distributionof the events recorded on the iteration (t):

$\begin{matrix}{{{f_{G}^{(t)}\left( {x,t} \right)} = {\sum\limits_{k = 1}^{\kappa^{*}}{w_{k}{f_{}\left( x \middle| {\overset{\sim}{Z}}_{k} \right)}{f_{{\overset{\sim}{Q}}_{k}}(t)}}}},} & (6)\end{matrix}$

The space-time density corresponding to all the events, which are notnecessarily observable, is obtained by normalization.Similarly, the entire distribution of the functional volumes isaccessible. We marginalize over time and define the spatial distributionof a functional volume j associated with the kinetic Q*_(j) on theiteration (1) by using the expression:

$\begin{matrix}{{f_{{VF}_{j}}^{(t)}(x)} = {\sum\limits_{{k:{\overset{\sim}{Q}}_{k}} = Q_{j}^{*}}^{\kappa^{*}}{w_{k}{{f_{}\left( x \middle| {\overset{\sim}{Z}}_{k} \right)}.}}}} & (7)\end{matrix}$

All the elements are then available to determine the estimators ofspace-time activity distribution and functional volumes VF for Titerations of the algorithm that are retained. The following expressionscan then be used for that:

$\begin{matrix}{{{\left( {\left. {f_{G}\left( {x,t} \right)} \middle| Y \right.,T} \right)} \approx {\frac{1}{T}{\sum\limits_{t = 1}^{T}{f_{G}^{(t)}\; \left( {x,t} \right)}}}}{and}} & (8) \\{{\left( {\left. {f_{{VF}_{j}}(x)} \middle| Y \right.,T} \right)} \approx {\frac{1}{T}{\sum\limits_{t = 1}^{T}{f_{{VF}_{j}}^{(t)}\; (x)}}}} & (9)\end{matrix}$

Similarly, any other functional of the a posteriori distributions can becomputed, like a conditional variance or a credibility interval.

1. A method for determining functional volumes for the search forkinetics representing the trend of the concentration of a radioactivetracer in an area of biological tissue, the method being applied tospatial components and comprising the following steps appliediteratively according to a Markov chain Monte-Carlo scheme: a step ofgeneration of a set K₀ ^(λ) made up of a set of candidate kineticsassociated with values of probability of appearance of these kinetics,these values depending on the concentration λ of radioactive tracer; alabeling step during which, for each spatial component of index K, theprobabilities of selection of the kinetics of the set K₀ ^(λ) areweighted by introducing therein a function λ_(k) representative of theconcentration of radioactive tracer in this component so as to obtain aset of indicative values, an indicative value D_(k) designating thekinetic with which the spatial component K is associated; a step ofconstruction of functional volumes, a functional volume VF_(j) beingmade up of the set of the spatial components which share the sameindicative value D_(k).
 2. The method as claimed in claim 1, furthercomprising a step of computation of the averages of the functionalvolumes VF_(j) obtained during the iterative process.
 3. The method asclaimed in claim 1, comprising a step of determination of statisticalestimators concerning the functional volumes VF_(j) obtained during theiterative process so as to quantify the estimation uncertainty.
 4. Themethod as claimed in claim 3, in which the variance of the functionalvolumes VF_(j) obtained during the iterative process is determined. 5.The method as claimed in claim 3, in which a credible interval isdetermined.
 6. The method as claimed in claim 1, in which the set K₀^(λ) is obtained by the determination (100) of two sets of scalarvariables V_(j) and Γ_(j), these variables relating to kineticcomponents of index j, in which V_(j) represents a fixed weightingcoefficient independent of the concentration λ and Γ_(j) represents areference concentration value associated with the kinetic of index j. 7.The method as claimed in claim 6, in which V_(j) is determined by usingthe expression:$\left. {a.\mspace{14mu} V_{j}} \middle| D \right.,{T\overset{ind}{\sim}{{Beta}\left( {{1 + m_{j}},{\nu + {\sum\limits_{l = {j + 1}}^{J_{\kappa_{n}}}m_{l}}}} \right)}}$an expression in which: for any j, m_(j) designates the number ofcomponents of the distribution H allocated to the j^(th) component of K₀^(Λ), H being a random measurement distributed according to a partiallydependent Dirichlet process, K₀ ^(Λ) being an innumerable collection ofrandom probability measurements which accepts, as distribution, a localDirichlet process; J_(k) _(n) designates the greatest index of thecomponents of K₀ ^(Λ) having at least one component of H allocated;Beta(,) represents the beta probability law.
 8. The method as claimed inclaim 6, in which Γ_(j) is determined using the expression:$\left. \Gamma_{j} \middle| V_{j} \right.,D,w,{\overset{\sim}{Z}\overset{ind}{\sim}{_{j}^{*}\left( \Gamma_{j} \right)}}$with${_{j}^{*}\left( \Gamma_{j} \right)} \propto {\left( {a_{j,\Gamma}^{\prime} < \Gamma_{j} < b_{j,\Gamma}^{\prime}} \right) \times {\prod\limits_{k \in _{j}^{+}}\; \left( {1 - {V_{j}\left( {{d\left( {\lambda_{k},\Gamma_{j}} \right)} < \psi} \right)}} \right)}}$expressions in which: Γ_(j) designates the parameter characteristic ofthe concentration associated with the kinetic j; D designates the set ofvariables D_(k) which indicate with which kinetic of K₀ ^(λ) the spatialcomponent K is associated; D_(j)−{k≦k_(n): D_(k)−j} represents the setof the indices K of the spatial components associated with the kineticj; D_(j) ⁺={k≦k_(n): D_(k)>j} represents the set of the indices K of thespatial components associated with a kinetic of index greater than j; wdesignates the set of the weightings w_(k) of the spatial components kin the Dirichlet mixture H; {tilde over (Z)} designates the set of theparameters, average and covariance matrix, {tilde over (Z)}_(k)associated with the spatial component K of the Dirichlet mixture H;

here indicates an independent random generation for each parameterΓ_(j); ∝ here signifies “proportional to”;

( ) represents the “gate” function such that

(condition)=1 if condition is true, and 0 otherwise; V_(j) designates acoefficient associated with each kinetic involved in the weighting ofsaid kinetic in K₀ ^(λ); d( ) represents the Euclidian distance in

: d(x₁,x₂)=|x₁−x₂| in which x₁, x₂∈

; ψ is a threshold chosen a priori which makes it possible toparameterize K₀ ^(λ).
 9. The method as claimed in claim 1, in which thelabeling step (109) is performed: by independently generating (500) forany k≦k_(n):$\left( {\left. D_{k} \middle| V \right.,v,Q^{*},\Gamma,w,\overset{\sim}{Z},C,T} \right)\overset{ind}{\sim}{\sum\limits_{j = 1}^{J^{*}}{p_{j,k}{\delta_{j}( \cdot )}}}$with${p_{j,k} \propto {\left( {{p_{j}\left( \lambda_{k} \right)} > \upsilon_{k}} \right){\max \left( {{p_{j}\left( \lambda_{k} \right)},\varsigma} \right)}{\prod\limits_{\{{1:{C_{i} - k}}\}}\; {{f_{Q_{i}^{*}}\left( T_{i} \right)}{and}{\sum\limits_{j = 1}^{J^{*}}p_{j,k}}}}}} = 1$by independently generating (501) for any K, such that k_(n)<k<k*:$\left( {\left. D_{k} \middle| V \right.,v,\Gamma,w,\overset{\sim}{Z}} \right)\overset{ind}{\sim}{\sum\limits_{j = 1}^{J^{*}}{p_{j,k}{\delta_{j}( \cdot )}}}$with p_(j, k) ∝ (p_(j)(λ_(k)) > υ_(k))max (p_(j)(λ_(k)), 𝜍) and${\sum\limits_{j = 1}^{J^{*}}p_{j,k}} = 1.$ in these expressions, k_(n)designates the number of spatial components with which at least onecoincidence event is associated; k* designates the total number ofspatial components of the Dirichlet mixture H for a given iteration; Vdesignates the set of the coefficients V_(j) involved in the weightingsof the kinetics mixture K₀ ^(λ); v represents the set of the auxiliaryvariables v_(k) used in the representation of the Dirichlet kineticsprocess; Q* designates the set of the Pólya trees Q_(j)* characterizingthe kinetics; Γ designates the set of the parameters Γ_(j)characteristic of the concentration associated with the kinetic j; Cdesignates the set of the classification variables C_(i) which associateeach coincidence event i with a spatial component of the Dirichletmixture H; T represents the set of the times of occurrences T_(i) of thecoincidence events; p_(j,k) designates the probability of associatingthe component K of the Dirichlet mixture H with the kinetic j; p_(j)(λ)represents the weight function taken at λ associated with theprobability of the kinetic j of K₀ ^(λ); δ_(j)( ) represents the Diracmeasurement function such that δ_(j)(j′)=1 if j=j′ and is 0 otherwise;

here indicates an independent random generation for each classificationvariable D_(k); max( ) represents the “maximum of” function; ζdesignates is a threshold parameter chosen arbitrarily involved in theimplementation of the slice sampling algorithm of the local Dirichletkinetics process; {i: C_(i)=k} designates the set of the coincidenceevents i which are associated with the spatial component K; f_(Q*) _(j)represents the kinetic j, i.e. the density of the Pólya tree j of thelocal Dirichlet mixture of kinetics; J* designates the total number ofcomponents of the local Dirichlet mixture of kinetics for a giveniteration.
 10. The method as claimed in claim 1, in which theconcentration λ_(k) is obtained by the application of a function φ_(λ)defined such that λ_(k)=φ_(λ)({tilde over (Z)}_(k), S) with:$\lambda_{k} = {{{_{{\overset{\sim}{x}}_{k}}\left( {\sum\limits_{m = 1}^{\infty}{w_{m}{\left( {\overset{\sim}{x}}_{k} \middle| {\overset{\sim}{Z}}_{m} \right)}}} \right)}\mspace{14mu} {with}\mspace{14mu} {\overset{\sim}{x}}_{k}} \sim {\left( {\overset{\sim}{x}}_{k} \middle| {\overset{\sim}{Z}}_{k} \right)}}$in which:

_(k) _(n) designates the mathematical expectation relative to theGaussian probability law of {tilde over (x)}_(k), N({tilde over(x)}_(k)|{tilde over (Z)}_(k)); N( ) represents the Gaussian probabilitylaw.
 11. A positron emission tomography device implementing the methodas claimed in claim 1.